clc
clear

cd 'D:\data_replication'

load('estimation/5_supply_side/main_data_supply.mat')
clearvars -except index_price_coeff

load('data/number_firms/companies_k.mat')

for i_k = 1:4092
    
i_k   
k = index_price_coeff(i_k,1);   

companies_k = company_summary{k};
N = table2array(companies_k(:, 4));
N(N < 1) = 1;
N_d = max(N);
N_d = round(N_d);
N_d = N_d * 2 + 100;

gamma_i_k = readtable('data/firm_size_distribution/gamma_i_k.csv');
gamma_i_k = table2array(gamma_i_k(:,2));

zeta_i_k = readtable('data/firm_size_distribution/zeta_i_k.csv');
zeta_i_k = table2array(zeta_i_k(:,2));

zeta_i_k(zeta_i_k < 0.7) = 0.8344;

prod_pct = ones(N_d, 1); 
for ii = 2:N_d
d_ii = log(ii) - log(ii - 1);    
%d_ii = (ii - (ii - 1))/(ii - 1);   
%d_prod_pct = - zeta_i_k(i_k,1) / 100 * d_ii;
d_prod_pct = - (1 / zeta_i_k(i_k,1)) * d_ii;
prod_pct(ii, 1) = prod_pct(ii-1, 1) * (1 + d_prod_pct);
end

prod_pct_summary{i_k} = prod_pct;


k = gamma_i_k(i_k,1);
a_min = 1;
n = 150000;
rng(3)
x = rand(n,1);
%a = k*(1-x).^(-1/a_min);
a = a_min*(1-x).^(-1/k);
a = sort(a, 'descend');
a = a(1:N_d,1);
prod_a = a / a(1,1);

prod_a_summary{i_k} = prod_a;

clearvars -except i_k prod_pct_summary prod_a_summary index_price_coeff company_summary

end

save('data/firm_size_distribution/prod_summary.mat', 'prod_pct_summary', 'prod_a_summary')


